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Abstract. We study numerically the "analyticity breakdown" tian- 
sition in l-dimensional quasi-periodic media. This transition cor- 
responds physically to the transition between pinned down and 
sliding ground states. Mathematically, it corresponds to the solu- 
tions of a functional equation losing their analyticity properties. 

We implemented some recent numerical algorithms that are ef- 
ficient and backed up by rigorous results so that we can compute 
with confidence even close to the breakdown. 

We have uncovered several phenomena that we believe deserve 
a theoretical explanation: A) The transition happens in a smooth 
surface. B) There are scaling relations near breakdown. C) The 
scaling near breakdown is very anisotropic. Derivatives in different 
directions blow up at different rates. 

Similar phenomena seem to happen in other KAM problems. 

Quasi-periodic solutions, quasi-crystals, hull functions, KAM theory 
[2000] 70K43, 52C23, 37A60, 37J40, 82B20 

1. Introduction 

In this paper we present detailed numerical results on the analytic- 
ity breakdown transition in models of one-dimensional quasi-periodic 
media. This analyticity breakdown transition is widely accepted to 
correspond to the transition between pinned down and sliding ground 
states. Hence, the position of the transition affects the properties of 
the material described by the models. For the physical motivation of 
the theory see Section |2| and |AA80llvE99] . 

From the mathematical point of view, the problem is to study the 
analyticity properties of the solution of a functional equation depend- 
ing on parameters (the function is the hull function parameterizing a 
ground state and the functional equation is the expression that the sys- 
tem is in equilibrium). It is expected that, for some parameter values, 
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the solution is analytic but that, as we move parameters, the domain 
of analyticity decreases and eventually disappears. See Section |2] for a 
more precise description of the problem. 

The phenomenon of analyticity breakdown is very similar to the 
widely studied phenomena of breakdown of KAM tori in mechanics. 
Indeed, the problem of analyticity breakdown in periodic media is com- 
pletely equivalent to the breakdown of KAM circles in two dimensional 
twist mappings. In the case of quasi-periodic media considered here, 
the problem does not seem to have an easy dynamical representation 
and the treatment of the functional equations describing the equilib- 
ria for KAM theory is very different from the Hamiltonian counterpart 
because, besides the lack of a dynamical interpretation, the system has 
an extra frequency. The KAM theory for the equilibrium equation in 
quasi-periodic media has been obtained recently in |iSdlL12at ISdlL12b] . 

The KAM theory in |SdlL12a] ISdlL12b] is very different from the 
regular KAM theory for two reasons: one is that the models do not 
admit an easy dynamical interpretation and a second reason is be- 
cause the quasi-periodicity of the substratum introduces the an extra 
variable. This extra frequency has very drastic effects, which are not 
just mathematical artifacts. For instance, there are counterexamples 
to some standard predictions of the Aubry-Mather theory which de- 
scribes the ground states in periodic media |LS03j . In the periodic 
case, when Aubry-Mather theory applies, the analyticity breakdown 
transition could be described as saying that the smooth tori break up 
into Cantori (see |Per79j ). In the case of quasi-periodic media, we are 
not sure of what could take the place of the Aubry-Mather theory after 
the breakdown of KAM theory. This is a question that deserves further 
study, but which we are not considering in this paper. 

Nevertheless, we expect that some of the phenomena studied in detail 
here, notably the anisotropic regularity at breakdown happen also in 
other KAM problems involving two-dimensional functions. Indications 
to that effect have already been found qualitatively |Tom96[ IHS99[ 
ICFL04j . but in this paper we can obtain rather quantitative aspects of 
the phenomenon. 

Since the main goal is to compute accurately functions as close as 
possible to their breakdown, it is clear that the numerics are going to 
be delicate and that it is very important to have criteria that allows to 
be confident that the computed solutions are correct. We note that the 
experience shows that close to the breakdown, the truncated equations 
admit many spurious solutions that do not correspond to truncations 
of true solutions of the equations |vEF02l IGvEOSj . 
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In this paper we have implemented the algorithms suggested in |SdlL12aj . 
The main results of |SdlL12a] are based on a Newton method which is 
very efficient (see Section 3.4). If the function is discretized using 
Fourier coefficients, the Algorithm |3. 2| gives quadratic convergence but 
requires only 0{N) storage and 0{N log{N)) arithmetic operations. 
Note that, even if we obtain the Newton-like quadratic convergence 
a step does not require to store (and much less to invert) an x 
matrix. 

Furthermore, the main result of |SdlL12a] is a result in a posteriori 
format. It asserts that, if we have a function which solves approxi- 
mately the invariance equation (e.g. a numerical solution of the equa- 
tion, where the error of truncation and round off is small) and that 
satisfy some easy to verify non-degeneracy assumptions corresponds 
to a true solution of the problem. Hence, by computing with a trun- 
cation, estimating the truncation error by changing the degree of the 
truncation and monitoring the non- degeneracy conditions, we may be 
reasonably confident that we have obtained a true solution. 

Of course, a numerical implementation requires many more details 
than a mathematical algorithm and one has to discuss implementation 
issues such as data structures used, choice of iteration steps, order of 
truncation, etc. These considerations are undertaken in Section |4j The 
results are presented in Section |4.1| and the conclusions are presented 
in Section [5] We anticipate that the conclusions are A) The breakdown 
happens in a smooth surface, B) There are scaling relations. Sobolev 
norms of the function blow up as a power of the distance to the critical 
surface. C) The breakdown is extremely anisotropic. 

Of course. A) and B) have been observed many times in phase transi- 
tions and are one of the predictions of renormalization group theory but 
C) seems somewhat unexpected. There have been visual observations 
indicating that the tori near breakdown look "filamented" . We think 
that there are reasons to expect that this phenomena will happen in 
other KAM problems involving several variables and we hope to come 
to that problem. In this paper, we will just present some numerical 
observations in a class of models. 



2. The problem of analyticity breakdown 

The simplest way to motivate the models studied here is to consider 
the problem of deposition of a material over a one- dimensional quasi- 
periodic substratum. 

If Ui represents the position of the i^^ particle of the deposited mate- 
rial, the state of the system is given by the configuration {ui}i^z, which 
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is sequence giving the position of all the deposited particles. Note that 
Ui will be a real variable. 

The Physics of the problem is obtained by assigning an formal energy 
to the configurations |Rue69t IIsr79t IMatOQj . This formal energy is 
obtained as a formal sum of the energies associated to every finite set 
of particles. Hopefully, most of the terms will be zero or decrease very 
fast with the diameter. Nevertheless, due to the translation invariance 
one does not hope that the sum of the energies converges. 

The main example considered in this paper will be the quasi-periodic 
Frenkel-Kontorova model given by: 

(2.1) S{{u}i(.i) = ^ ^{Un - Un+lf - V{aUn) 

where : T'^ — )■ M and a G M'^ is a sufficiently irrational vector. 

The term \{un — Wn+i)^ represents the interaction between nearest 
neighbors and the term V{aUn) represents the interaction with the sub- 
stratum at the position «„. Note that this function is quasi-periodic as 
a reflection of the fact that the substratum is quasi-periodic. The clas- 
sical Frenkel-Kontorova model jFK39] considered the case of periodic 
potentials with d = 1. 

Other physical interpretations of the model are possible, for exam- 
ple, the classical |FK39] introduced the Frenkel-Kontorova model to 
describe the motion of planar dislocations in a crystal. 

2.1. Some standard definitions. Now we recall some standard def- 
initions |AA80lRLD83] . 

A configuration {Mijiez is in equilibrium when 

(2.2) du,S{{ui},^j) = 



In the case (2.1), the equilibrium equations (2.2) become: 



(2.3) a ■ VV{aUi) + Wj+i -|- — 2ui = 



Note that, even if the sum in (2.1) is a formal sum, the equilibrium 



equations (2.3) are very well defined. 



A configuration {ui}i^z is a ground state when 
(2.4) S{u + 7]) - S{u) >0 

for all 7] that vanish at all but a finite number of indices. Note that the 



left hand side of (2.4) can be given a good interpretation in the model 
(2.1) because the two formal sums only differ in a a finite number of 
sites. 
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A configuration is given by a hull function when we can find h : 
T'^ ^ M and w G M such that 

(2.5) Un = (jjn + h{nuja). 

2.2. Equilibrium equations for hull functions. Formulation of 



our problem. Substituting (2.5) into (2.2), we obtain that an equi- 
librium configuration is given by a hull function if and only if 
(2.6) 

h{nuja + ua) + h{nuja — uoa) — 2h{nu)a) + (a ■ V) V{nuja + h{nuja)). 

If {nua}n(^z is dense on T'^ (which is well known to happen if and 
only if k ■ ua ^ Z for /c G Z'^ — {0}, [KH95j), for a continuous h we 



obtain that (2.6) is equivalent to 

(2.7) h{e + ioa) + h{e - ooa) - 2h{e) + (a ■ V) ¥{6 + h{e)). 

Of course, it could well happen that there are discontinuous hull 
functions, but we are interested in studying precisely the existence of 
continuous - indeed analytic - functions. 



Hence, (2.7) will be the centerpiece of our attention. We will fix oja 
satisfying good number theoretical properties, assume that V depends 
on parameters and study numerically the set of parameters for which 



(2.7) has an analytic solution. We will pay special attention to the 



behavior of these solutions near the boundary of existence. 



We note that the existence of a continuous solution of (2.7) implies 
that there is a continuous family of equilibrium configurations. On can 
note that for any ao eT'^, ha-^ given by 

(2.8) h,^{0) = h{e + cxo) 



is also a solution of (2.7). Hence, the configuration corresponding to it 
according to (2.5) is also an equilibrium configuration. The fact that 
we obtain a continuous family of configurations which are equilibria 
indicates that there is no energetic impediment to the solutions making 
small transitions and sliding. Hence one expects that the solutions can 
be modified easily. In contrast, if the set of equilibrium configurations 
is discontinuous, there may be an energetic barrier (the Peierls-Nabarro 
barrier) to perform jumps from a configuration to the nearest one and 
the solutions are pinned down. 

Notice that hcrg corresponds to choosing a different origin of coor- 
dinates to the internal phase of the problem. This can be considered 
as a gauge symmetry of the problem. The Ward identities associ- 
ated to this gauge symmetry play an important role in the study in 
jdlLOSl ICdlLlObl ISdlL12al ISdlL12bj and also in our numerical treat- 
ment. 
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3. Mathematical results 
3.1. An overview of the rigorous results. In the periodic case (rf 



1), the problem of existence of continuous hull functions satisfying (2.7) 
is equivalent to the problem of existence of rotational invariant circles 
for a twist mapping |ALD 83]. In the periodic case, there are several well 
developed mathematical theories which can lead to an understanding 
of the problem. 

• Kolmogorov- Arnold-Moser (KAM) theory |Mos66llHer83] . 

• Aubry-Mather theory |ALD83l [Mat82l IMos86j . 

• Renormalization group theory |Mac93t IKoc04j . 

In the periodic case, there are a number of numerical methods to study 
the breakdown. The appendix A of [CdlLlOb] contains a comparative 
summary of the literature and compares the methods. 

The KAM theory is perturbative but produces analytic solutions of 



(2.7) when the frequency satisfies some number theoretic properties 
(called "Diophantine") . Aubry-Mather theory is non-perturbative and 
it does not require any number theoretic properties of the function. 



On the other hand, the solutions of (2.7) may be discontinuous. The 
renormalization theory has as an aim to describe a subset of the bound- 
ary of analyticity. It was originally developed based on non-rigorous 
arguments supported by very careful numerics |Mac93] . but by now 
there are quite a number of rigorous results (some of them proved by 
computer-assisted proofs) . 

In the quasi-periodic case considered here, the situation is very differ- 
ent. The Aubry-Mather theory seems somewhat problematic due to the 
examples in |LS03] (but there are some topological results in |GGP 06]). 
To the best of our knowledge, there is no renormalization group theory 
for general quasi-periodic FK models (see, however |MQ00j ). There 
has been a recent development of a KAM theory in |SdlL12a] ISdlL12b] 
which will be the basis of our work. 



3.2. A posteriori results. The main result of |SdlL12a] is formulated 
in an a-posteriori format standard in numerical analysis. 
We define the operator T acting on hull functions by: 

(3.1) r[h]{e) = h{e + ua) + h{e - ua) - 2h{e) + (« ■ v) v{e + h{e)) 



So that (2.7), the equilibrium equation for hull functions can be con- 



cisely expressed as 



(3.2) 



r[h] = 
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The main result of |SdlL12a] says (omitting some mathematical as- 
sumptions on regularity, definition of norms and number theoretic prop- 
erties of Loa) that if we can find a function that satisfies approximately 
the equation (in some norm) and satisfies some non- degeneracy condi- 
tions, we can find a true solution and bound the distance from the true 
solution to the original guess (measured in another norm) by a constant 
times the size of the residual of the initial approximate solution. 

Prototype Theorem 3.1. Assume that we can find a function : 
^ M. such that 

lirNlli<e 

Mi(/io) < Ai,...,M„(/io) < A 



m 



then, if 

e < e*(Ai,...,v4^) 

there exists a true solution h* ofT{h*) = 0. and 

\\h* - ho\\2< K{Ai,--- ,Am)e 
Furthermore, the solution is the only solution (up to the gauge trans- 



formations associated to the change of origin in (2.8) ) 



In the above. Mi, . . . , are rather explicit condition numbers, 
and e* is also an explicit function. Note that, as customary in KAM 
theory, there may be a loss of regularity and that the norm in which we 
measure the error may be different than the norm in which we reach 
the conclusions. In jSdlL12a] . one can find results where the norms 
II ■ 111, II ■ II2 are norms in spaces of analytic functions and also results 
for Sobolev norms. 

In applications, the approximate solution will be the product of a 
numerical calculation. It is important that in view of the above result, 
that to be confident that the approximate solution corresponds to a 
true solution, we do not need to study the algorithm used to produce 
it. We just need to verify that the condition numbers are reasonable 
and that the numerical error is small compared to them. 

As we will see, the method of proof in |SdlL12aj is based on a rapidly 
convergent iterative method which is explicitly described. This method 



leads to an efficient algorithm (See Section 3.4, which we have imple- 
mented in Section |4] and which is the basis of our results. It is also 
important that the method leads to a numerically accesible criteria for 



breakdown, which we discuss in Section 3.3 
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3.3. A numerically accessible criterion for breakdown of an- 
alyticity. In |(MlLinbj . it was shown in great generality that, when 



one has an a posteriori result of the form Prototype 3J^ working for 
Sobolev and analytic spaces, one can get several results automatically. 

Bootstrap of regularity All the solutions which are in a Sobolev 
space of order r > tq are analytic. 

The critical value tq depends only on the number theoretic properties 
of ua. 

Criterion for breakdown of analyticity If we consider a family 
depending on parameters, and obtain a family of approximate solutions 
with bounded non-degeneracy conditions, the parameters approach the 
boundary if and only if Sobolev norms of high enough order go to 
infinity. 

Hence, it is clear how one should proceed. Implement the algorithm 



described in Section 3.4, run it monitoring the condition numbers (so 
that we are sure that we are not considering any spurious solutions) 
and the Sobolev norms. If the Sobolev norms blow up dramatically 
and nothing else happens, we are approaching the breakdown. Since 
the algorithm gives the Fourier coefficients of the hull function, the 
computation of the Sobolev norm is rather straightforward. 

It is important to note that this criterion does not require a dynam- 
ical interpretation and in |CdlL09l ICdlLlOb] one can find a study of 
systems with long-range interaction for which there is no dynamical 
interpretation. The paper [CdlLlOb] contains a comparison with other 
methods. 

The above criterion for blow up has been implemented several times. 
For twist mappings it has been implemented in |CdlL09[ ICdlLlOa 



IFM12j . For dissipative systems it has been implemented in |CC10t 
[CFT2] . 

3.4. An efficient algorithm for a Newton method. The proof of 
|SdlL12a] is based on a Newton-like method that given an approximate 
solution produces a much more accurate one. The error after the cor- 
rection is, roughly, the square of the error before the correction, but 
there is a loss of regularity. It is well known in the Nash-Moser hard 
implicit function theorems that these methods still lead to convergence. 
The algorithm in [SdlL12a] , which we detail below, is obtained taking 



advantage of several identities satisfied by (2.7) that come from its 
variational character and the gauge symmetry (2.8). 

These identities reduce the Newton step to a sequence of elementary 
operations on functions (composing, taking derivatives, multiplying, 
solving difference equations with constant coefficients). A remarkable 
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feature is that all of these operations are diagonal either in a discretiza- 
tion in Fourier terms or in a discretization based on points on a grid. Of 
course, the Fast Fourier Transform allows to pass from one discretiza- 
tion to another in a easy way. 

The upshot is that if we discretize the function h in terms of N 
Fourier coefficients and the values at N points in a grid, a Newton- like 
step requires 0{N) storage and 0{N\og[N) operations. 

In practice, the constants in front of the asymptotics are moderate, 
since there are very optimized implementations of FFT for almost any 
architecture. As we will see, it is quite possible to use N ^ 10^ in a 
common desktop. 

For the sake of completeness, we include verbatim the algorithm from 
|SdlL12a] . 

The algorithm of jSdlL12a] is actually slightly more general instead 



of (3.2) it consideres the equation 

(3.3) T[h] + A = 

where A G M, and both h and A are the unknowns. 



The equation (3.3) allows to consider also also forces that do not 



come from a potential. Nevertheless, it is shown in |SdlL12a] that 



when the forces come from a potential, any solution of (|3.3|) satisfies 
A = and, hence gives a solution of ( 3.2[ ). 



As we will see, the algorithm consists of steps which are all elemen- 
tary manipulations of functions (derivatives, multiplications, compo- 
sitions) as well as solving the so called cohomology equations. That 
is, given h{9) periodic of period 1 and with j h{9) dO = 0, find W also 
periodic with period 1 and with average such that 

(3.4) W -W oT^ = b 

where T^(^) = e + a. 



Note that (3.4) can be solved very efficietly using Fourier coefficients 



and that it is equivalent to the equation for Fourier coefficients 

The factor in parenthesis in the LHS vanishes for k = 0, and the 
equation for = is overdetermined. It does not have a solution if 
bo 7^ 0, but if 6o = any Wq is a solution, we choose the solution with 
Wo = 0. If w is irrational, the term in parenthesis does not vanish 
for k ^ and, if it is Diophantine, we have lower bounds for it that 
lead to estimates for the solutions of the cohomology equation |Riis75t 
IRiis76j . These estimates are used in |SdlL12aj to prove convergence of 
the algorithm, but in this paper, we will not discuss estimates. 
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Algorithm 3.2. Given h:T'^ ^R, A G M with h{e) = Y.k&i he^^'''^ 
and h(t) = h{at) for t G M and any irrational vector a G M'^, we will 
perform the following calculations (where (■) denotes the average) 

1) Let ^ = h{t + uj) + h{t -uj)- 2h{t). 

In Fourier components =Sffc = 2(cosco'a ■ k — l)hk- 

2) We calculate U = V{e + h{e). 

3) Calculate e = ^ + U + \ = T[h] 

4) Calculate I = 1 + dah. 

In Fourier components Ik = Sk.o + 2TTi{k ■ a)hk- 

5) Let f = i-e. 

6) Choose 5 = — /q. 

7) Denote b = i-{e + S). 



8) Solve the cohomology equation (3.4) for W with zero average. 



That IS, Wj, - 2(cos..».fc-i) - 
9) Take W = 



10) Calculate W = W° + W . 

11) Solve the cohomology equation for A . That is, = 2(cos "q fc-i) 



where a 



w 



12) We obtain A = A ■ /. 

13) The new improved solution is h + A, X + S. 



4. Numerical Implementation 



The most computationaUy expensive steps in Algorithm |3.2| are di- 
agonal in either real space or in frequency space, so we use the Fast 
Fourier Transform multiple times in each application of the Newton 
step. We use Matlab, and typically 2^^ Fourier modes, so a discretiza- 
tion for his a 512 x 512 complex-double. All non-degeneracy conditions 
are satisfied if the norm of the error is small enough, so we monitor the 
error throughout. 

To implement the algorithm, we compute the error of an initial guess, 
then perform the Newton step to improve the guess, and then compute 
the new error for the improvement. If the new error is not below a 
threshold, the Newton step is performed again and the process repeats. 
If the new error is below the threshold, we conclude a torus exists 
nearby and we increase the parameter values. If the Newton step fails 
to bring the error below the threshold (after some number of attempts) 
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or if the error increases, we end the program and start again with a 
better initial guess. 

Wc explore the parameter space along rays at a fixed angle 9 by 
starting at {61,62) = (0,0), and increasing {61,62) — )■ (^i + Ae, £2 + 
A£tan6'). We use the most recently computed invariant torus as our 
initial guess (possibly scaled a bit) at the next parameter value. If the 
Newton step cannot bring the error low enough after a certain number 
of iterations, we take data from the last point of convergence, and try 
again but with a close parameter value. For example, if the program 
converged with (£1,^2) and we take the step 61 ^ 61 + A^i, 62 
£2 + A£2, at which point the program fails, then we try again with 
Si ^ 61 + ^Asi, 62 ^ 62 + |A£:2- If the program fails to converge, and 
we have taken a very small step in parameter space, then we terminate 
the program and consider the lack of convergence as an indication that 
the boundary of analyticity is very close. 

The Newton step itself can bring in errors. The FFT is used exten- 
sively, and occasionally a mode hk with small modulus will appear from 
computer error, which can be amplified as the iterations continue. This 
can lead to problems with convergence due to artificially large Sobolev 
norms. We have found it necessary to monitor the appearance of false 
modes and to eliminated them. We implement a method that sets to 
zero all modes whose amplitude is below a certain threshold: they are 
deemed to be numerical errors and are set to zero so no amplification 
is possible. This slows the convergence of the Newton method, but 
allows us to compute tori much closer to the boundary of analyticity. 
This is a delicate process; if done poorly, it can cause the Newton step 
to return the zero function, and thus fail to find an improvement. We 
use a cut-off threshold that is allowed to vary. 

When implementing the algorithm, we can only consider a finite 
number of Fourier modes, which can lead to "spurious" solutions. That 
is, solutions of the discretized problem that do not correspond to so- 
lutions of the original problem. For example, when we are looking for 
a solution to T[h] = 0, and we replace this with a discretized problem 
Thlh] = using N Fourier modes, we need to ensure the discretized 
solution /ijv satisfying T/v[^'iv] = will also satisfy Tlh^] ^ 0. To avoid 
computing discrete solutions that satisfy 7iv[/iAr] = but do not 
satisfy T[/iAr] ~ 0, we compute the error of on a finer grid. That is, 
check if 72iv[^Ar] ~ 0. 
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We will consider two different models, whose potentials are given by 
(4.1) 

^1 </'2) = cos(27r</)i) - — cos(27r(/)2), 
(4.2) 

^2(01, 02) = cos(47r0i + 47r02) - tt- (cos(27r0i) + cos(27r02)) • 
47r zvr 



4.1. Results for Models (4.1) and (4.2). The results we present are 
for the values 



(4.3) 



a 



1.246979603717467, 2.801937735804838). 



These numbers are roots of cubic polynomials 



X 



4x2 + 3x + l, 



a; — 2x — 1), but we found similar results for other numbers. 

For a fixed rotation, we use the algorithm to explore rays in param- 
eter space {ei,e2)- In Figure [l] we plot the domains of analyticity for 
Models ( 4.1[ ) and (4.2), that is, values of {61,62) where the algorithm 
converged. The boundary of the domains of analyticity for each model 
appear to be smooth curves (at least locally). For Model (4.1) the 



domain appears to have reflection symmetries over both the and £2 
axes. Model (4.2) appears to only have reflection symmetry over the 
El axis. 



Figure |2] shows an invariant torus (from Model (4.1)) for parameter 
values near the breakdown of analyticity (2a and 2b are two different 



views of this torus). One can see that the torus oscillates rapidly in a 
single direction (this is the direction ua). This "filamented" appear- 
ance indicates an anisotropic breakdown: the derivatives of h blow up 
faster in the direction ua than they do in the direction ua-^. Figure |3] 
shows two curves that lie on the surface depicted in Figure |2j These 
curves are obtained by cutting the surface in Figure [2] along the di- 
rections ua and uja-^. Along, ua, we see fast oscillations (left graph) 
compared to the cut along ua-^ (right graph). (A more quantitative 
description of the anisotropic breakdown is given in Figure [5} described 
below.) 

In Figure |4] we have depicted the indices k for which the correspond- 



is far from zero. We 
brm is clustered along a 



ing Fourier mode hk of the torus in Figure 
observe that the support of the Fourier trans 
line in Fourier space. 

That means that the function h can be aproximated by a function 



A(9i,«2 



27rjA;6»i+/3fc6»2 



H{d, + 1392) 
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(a) Domain of analyticity for Model (4.1). 



X 10 




x 10 



(b) Domain of analyticity for Model (4.2). 



Figure 1. Values of {61,62) where the models have an- 
alytic invariant tori. 
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(a) Torus near breakdown for Model (4.1) 




(b) Rotated view of the torus. 
Figure 2. Two views of a torus near breakdown for 



Model (4.1). Ridges appear as though h has a frequency 



close to ua. 



which also lends credence to the belief that the result near breakdown 
is basically a one dimensional function. 

Notice that, once we know that the final result is essentially a one 
dimensional function, it should be possible to perform linear changes 
of the 6 variables so that this line is approximately parallel to one of 
the coordinate axis. Then, it would be more efficient to choose more 
Fourier modes along this direction. Of course, since our goal was to 
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(x, a /a x) 



Figure 3. Left plot is of a cut along the direction loa 
of the torus from Fig. |2} The right plot is a surface cut 
along the direction loa^ . The torus is much smoother 
along the uja^ direction. 



illustrate the appearence of this linear structure, we decided to use 
general algorithms that do not take advantage of its presence. 

In Figure ISj we plot a log-log scale of the ifij-norms and if^-norms 
(for r = 4, 5 j versus ei^crit — where ei,crit is the critical value of e\ 
along the line £2 = tan (f ) ^i- We see that the ifpnorms blow up much 
faster than the if^-norms, exhibiting the anisotropic breakdown. 

In Figure 6a, we plot ||/i||r,|| against ||/i||8,||, and ||/?.||r,± against ||/?.||8,_l 
for r = 9, 10 on log-log scales. As the parameters (^i, £2) approach crit- 
ical values, the log-log scale behavior is linear. (The if -norm behaves 
like the ifpnorm because of the anisotropic behavior.) 

The blow-up behavior seen in Figure [sjfits a power-law form ||/i||r,|| ~ 
C{ecx\\. — e)^*-''-', with exponent p depending on the order of the Sobolev 
norm. Computing the exponent for different values of r, we find a 
nearly linear dependence on r. In Figure 6b we plot the blow-up ex- 
ponent for values r = 1, 1.5, 2, . . . , 10 as computed along two different 
lines in parameter space for Model (4.1). The blow-up exponent is 
■p{r) = —f3\\r + 7|j for some constants (3\\, 7||. The fit-line in the left 
plot has (3\\ ?a 0.34 and 7|| ^ 0.41, while in right plot (3\\ ^ 0.33 and 
7ll ~ 0.26. With higher-accuracy numerics, we expect to have similar 
plots for the iY^-norms, too. These subdominant norms are measured 
with lower accuracy and more susceptible to noise in the numerics. We 
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250 



frequency 



Fourier modes for a torus before breakdown 
This is the same torus as in Figure 
A dot is plotted at the point k = {ki,k2) if Fourier 
coefficient hk is non-zero. Using 512^ coefficients, we have 
-256 < ki,k2 < 255. 



Figure 4. 
for Model ( [il] ) 
2 



can see that they blow up like a power as (£1,62) approach the bound- 
ary of analyticity (Figure |5]), but measuring the exponents as functions 
of the regularity will require more powerful numerics. 

5. Conclusions 

We have studied numerically the problem of breakdown of analyticity 
of hull functions in quasi-periodic media. This transition corresponds to 
the physical effect of transition between pinned down or sliding ground 
states. 

Our computation has been carried out by efficient algorithms backed 
up by a-posteriori theorems that allow to compute confidently ex- 
tremely close to the break up. We have uncovered certain regularities 
that, following the style of |Gre79] . we call Assertions. 

Assertion I The boundary of the domain of analyticity is a smooth 
curve (at least locally in the region explored). 

Note that, following |AA80 ] this is analogous to what happens in 
phase transitions, where the phase transitions - which also correspond 
to a breakdown of analyticity - happen often in smooth surfaces (Gibbs 
rule) |Isr79j . The fact that the transitions happen locally in a smooth 
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manifolds is one of the predictions of the renormahzation group because 
they are the stable manifolds of a fixed point. 

Of course, it can well happen that there are regions where the phase 
transition has cusps or other singularities, but we have not encoun- 
tered them. It is, of course, quite conceivable that a more detailed 
exploration will uncover these singularities. As a cautionary tail, we 
mention that the paper [CdlLlOb] uncovered singularities in the case 
of twist mappings that appeared only when a very large range of pa- 
rameters was explored. The singularities are an indication that the 
renormahzation group dynamics may have a more complicated behav- 
ior. Note that the apparent cusps found in Figure[T]near the coordinate 
axis are not real singularities. The potential has symmetries that the 
reflections along axis with the change in the origin. 

Assertion II There are scaling behaviors near the blow up 

We have seen that near the blow up we have 

where /3 (and perhaps 7 too) is a universal number. 

Actually, we have found a more precise version of scalings in each 
direction. 

Assertion III The blow up of the norms is very anisotropic 
If we define 

||/^||,,ll = ||(a;«-Vr/i||z.. 
\\h\\r,± = \\iua^-\/Yh\\L2 
we have near the breakdown 

Moreover, when e — Ec, we have: 

||^£||r,±/||^e||r,|| 0. 

Note that, as a particular case of the assertion, the tori remain some- 
what regular up to the breakdown along lines. 

Remark 5.1. We think that it is also possible that there are scaling 
relations 

with smaller exponents than those for the parallel norm. Unfortu- 
nately, verifying the above scaling relations is beyond reach at the mo- 
ment. The fact that the asymptotic expansion is subdominant makes it 
harder to compute numerically (the calculations are contaminated by 
the larger effect) and also one can expect that this effect will manifest 
itself only much closer to breakdown. 
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vs. £1 crit - ei for r = 4,5. 



Model (4.1). 
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(b) Log-log plots of ||/i||r.± VS. Gi^crit — El for r = 4, 5. Model 



Figure 5. ||/i||r,|| and ||/i||r,± blow up at different rates as 
61,62) increase along the line £2 = fan (|) ei for Model 



(4.1 ). 5a shows blow-up behavior for Model H-lh, and 5b 



shows blow-up behavior for Model (4.2). The ifiy-norms 



blow up much faster than the if^-norms. 



22 



T. BLASS AND R. DE LA LLAVE 



10 



10 



o 10 



10 



10 



10 



10' 



§ 10^° 



10 



10 



10'" 

logllhiL 



10'" 

logllhii, 



10 ' 



10 

10^" 10° 



10 



10 ' 



10 



10 10 



10'" 

logllhiL 



log I 



10 



/ 



10 



10 



(a) Log-log plots of 
r = 9, 10. 
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(b) Blow-up exponents for 
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VS. r. 



Figure 6. 



6a shows log-log plots of Hu vs Hu and vs 



^8 ^Qj^ r = 9, 10. (All plots in this figure come from tori 
for Model (4.1)). 6b shows the relationship between the 
blow-up exponent of \\h\\ri\ function of r as {81,62) 
move along the rays 62 = tan (|) £1 and £2 = tan (|) ei. 
We see that \\h\\r,\\ ^ c{e, - e)-'^ir+T|i. 



